Abstract
Resource planning for hospitals under special consideration of the COVID-19 pandemic.
A sector as critical as health sector experiencing unplanned demand for resources can lead to devastating loss. This challenge was experienced heavily around the world since the COVID19 pandamic started. Though, the challenge has always been there in the health sector, the pandemic amplified it.
Lack of suitable scientific tool and for modeling, predicting and planning for resources usage in the hospital sector has made the challenge more difficult to deal with. Hence, the need for a suitable scientic/data analysis tool to deal with this. This is the reason and idea of the Babsim Resource Planning tool.
The Babsim R package is tool for Simulation of resource allocation in hospitals.According to Lawton and McCooe paper on the same subject, “There are a variety of ways in which ICU capacity, and other hospital processes, can be modelled. These can be categorised into two main types: top-down (often ‘system dynamics’ or ‘time series’), and bottom-up (usually ‘discrete event’ or ‘agent-based’) approaches” The tool uses the second method, which is Discrete event approach to do the simulation.
A discrete-event simulation (DES) which is used in this simulator models the operation of a system as a (discrete) sequence of events in time. Each event occurs at a particular instant in time and marks a change of state in the system. Between consecutive events, no change in the system is assumed to occur; thus the simulation time can directly jump to the occurrence time of the next event, which is called next-event time progression.
SPOT: A set of tools for model-based optimization and tuning of algorithms.
Simmer: Simmer is a process-oriented and trajectory-based Discrete-Event Simulation (DES) package for R.
The required packages are installed if not already exists.As shown below, the packages are devtools for working with git, spot for working with spot and the babsim package.
#install.packages("devtools")
#devtools::install_github("r-lib/devtools")
#url <- "http://owos.gm.fh-koeln.de:8055/bartz/spot.git"
#devtools::install_git(url = url)
#url <- "http://owos.gm.fh-koeln.de:8055/bartz/babsim.hospital"
#devtools::install_git(url = url, subdir = "babsim.hospital")
rm(list = ls())
suppressPackageStartupMessages({
library("SPOT")
library("babsim.hospital")
library("simmer")
library("simmer.plot")
library("plotly")
})
The package version of SPOT is checked since not all version could be used in the program. Package version of SPOT must be larger than 2.2.4:
Package version of SPOT must be larger than 2.2.4:
packageVer
# Motivation
* Crises like the COVID-19 pandemic pose a serious challenge to health care institutions.
* They need to determine and plan the resources required for handling the increased load,
for instance in terms of hospital bedssion("SPOT")
and ventilators. After consultation with the local authorities of the Oberbergisches Land District, the babsimhospital tool was created to address the resource planning challenges the public authorities are facing. * This is a tool for capacity planning based on discrete event simulation that aims to address this challenge. * The predictive quality of the simulation is based on a set of 29 parameters. * The default values of these parameters were established in cooperation with medical professionals. * We aim to investigate and optimize these parameters to improve babsimhospital. * To that end, we use model-based optimization via SPOT and an in-depth sensitivity analysis. * The sensitivity analysis is crucial for the optimization process, since it allows to focus the search on the most important parameters of the simulation. * We illustrate that this approach allows for a reduction in the dimensions of the simulation without compromising the resulting accuracy.
We combine data from two different sources:
simData: simulation data, i.e., input data for the simulation. Here, we will use data from UK.fieldData: real data, i.e., data from the DIVI-Intensivregister. The field data will be used for validating the simulation output.Statistically speaking, the babsim.hospital simulator models resources usage in hospitals, e.g., number of ICU beds (\(y\)), as a function of the number of infected individuals (\(x\)).
In addition to the number of infections, information about age and gender can be used as simulation input.
You will take a closer look at these data in the following sections. Therefore, we have provided som real data from the UK.
library(readxl)
### X20201111UKdata <- read_excel("/Users/bartz/workspace/Lehre.d/IDEA-Master-AIT-WS2020-2021/Numerical-Methods/Datasets/CovidDataGroup1.xlsx")
X20201111UKdata <- X20201111UKdata <- read_excel("CovidDataGroup3.xlsx")
ukdataRawFull1 <- X20201111UKdata
str(ukdataRawFull1)
%> tibble [240 x 5] (S3: tbl_df/tbl/data.frame)
%> $ Date : POSIXct[1:240], format: "2020-03-16" "2020-03-17" ...
%> $ TotalCOVID19Inpatient : num [1:240] 0 2 3 0 3 10 13 19 22 34 ...
%> $ COVID19NonInvasiveCPAP: num [1:240] 0 0 0 0 0 0 0 0 0 0 ...
%> $ COVID19VentilatedICU : num [1:240] 0 0 0 0 0 0 0 0 0 0 ...
%> $ NewCases : num [1:240] 0 5 3 0 3 5 5 10 14 11 ...
Infected: New cases = new cases either admitted to hospital or diagnosed in the hospitalbed: Total COVID inpatient = total in hospital (includes ICU)ìntensiveBed: COVID-19 Non-invasive = patients on non-invasive ventilators (CPAP). In theory these should be on ICU but we don’t have space for them all. The model should probably consider them to be ICU but not intubated (level 2)intensiveBedVentilated: COVID-19 Ventilated = patients ventilated and intubated on ICUNewCasesUK: New cases in UK = new cases diagnosed across the city. In March/April the government testing was very poor so there were a lot of undiagnosed cases in the UK. Currently we think we’re diagnosing 30-50% of cases.ukdataRaw <- ukdataRawFull1[as.Date(ukdataRawFull1$Date) > "2020-09-01",]
str(ukdataRaw)
%> tibble [70 x 5] (S3: tbl_df/tbl/data.frame)
%> $ Date : POSIXct[1:70], format: "2020-09-02" "2020-09-03" ...
%> $ TotalCOVID19Inpatient : num [1:70] 18 21 21 19 22 19 18 19 21 22 ...
%> $ COVID19NonInvasiveCPAP: num [1:70] 3 5 5 5 5 5 6 5 5 5 ...
%> $ COVID19VentilatedICU : num [1:70] 2 2 2 2 2 2 2 2 2 0 ...
%> $ NewCases : num [1:70] 112 93 91 98 114 182 82 122 152 126 ...
simData <- data.frame(Day = as.Date(ukdataRaw$Date),
Infected = ukdataRaw$NewCases)
plot(simData$Infected ~ simData$Day, type="b")
Note: non ICU patients can be calculated as TotalCOVID19Inpatient - COVID19NonInvasiveCPAP - COVID19VentilatedICU
We convert the data into a data.frame with
Day <- as.Date(ukdataRaw$Date)
bed <- ukdataRaw$TotalCOVID19Inpatient - ukdataRaw$COVID19NonInvasiveCPAP - ukdataRaw$COVID19VentilatedICU
intensiveBed <- ukdataRaw$COVID19NonInvasiveCPAP
intensiveBedVentilation <- ukdataRaw$COVID19VentilatedICU
fieldData <- data.frame(Day = Day,
bed = bed,
intensiveBed = intensiveBed,
intensiveBedVentilation = intensiveBedVentilation)
str(fieldData)
%> 'data.frame': 70 obs. of 4 variables:
%> $ Day : Date, format: "2020-09-02" "2020-09-03" ...
%> $ bed : num 13 14 14 12 15 12 10 12 14 17 ...
%> $ intensiveBed : num 3 5 5 5 5 5 6 5 5 5 ...
%> $ intensiveBedVentilation: num 2 2 2 2 2 2 2 2 2 0 ...
max(fieldData$intensiveBed)
%> [1] 64
bed: non ICU patients in hospitalintensiveBed: ICU bed without ventilationintensiveBedVentilation: ICU bed with ventilationTo perform a simulation, the setting must be configured (seed, number of repeats, sequential or parallel evaluation, variable names, dates, etc.)
To perform a simulation, the setting must be configured (seed, number of repeats, sequential or parallel evaluation, variable names, dates, etc.)
Seed ensures the same random number is generated
parallel determines how the simulation will be run, parallel or not
percCores controls the utilization of the computer processor
the resourceNames is list of the names resources to be modeled
resourceEval is similar to the resourceNames, it shows the actual title of the resources to be modeled.
seed = 123
simrepeats = 2
parallel = TRUE
percCores = 0.8
resourceNames = c("bed", "intensiveBed", "intensiveBedVentilation")
resourceEval = c("bed", "intensiveBed", "intensiveBedVentilation")
Other configuration is performed, the start date for simulation is determined. The fielddata is initialized with null and the icu (requirement) is initialized as false. The Weight of different bed types is also specified. As can be seen, the resources weight of the intensive care bed is very high compared to others.
We can specify the field data based on ukdataRaw for the simulation as follows:
FieldStartDate = as.Date(min(fieldData$Day))
rownames(fieldData) <- NULL
icu = FALSE
icuWeights = c(1,2,10)
SimStartDate = FieldStartDate
data:data <- list(simData = simData,
fieldData = fieldData)
visualizeUK(data=data)
UKdata <- data
### usethis::use_data(UKdata, overwrite = TRUE)
Next we extract the simulation configuration using the babsimToolsConf function. As seen below, the configuration is the function can be changed. Later the configuration is stored in para variable and used in simmer.
Configuration information is stored in the conf list, i.e., conf refers to the simulation configuration, e.g., sequential or parallel evaluation, number of cores, resource names, log level, etc.
conf <- babsimToolsConf()
conf$ResourceEval <- conf$ResourceNames
conf <- getConfFromData(conf = conf,
simData = data$simData,
fieldData = data$fieldData)
conf$parallel = parallel
conf$simRepeats = simrepeats
conf$ICU = FALSE
conf$ResourceNames = resourceNames
conf$ResourceEval = resourceEval
conf$percCores = percCores
conf$logLevel = 0
conf$w2 = icuWeights
set.seed(conf$seed)
The simulations using simmer (Discrete Event simulation) package to simulate resource utilization in the hospital.
The core of the babsim.hospital simulations is based on the simmer package.
It uses simulation parameters, e.g., arrival times, durations, and transition probabilities.
These are currently 42 parameters (shown below) that are stored in the list para.
para <- babsimHospitalPara()
str(para)
%> List of 29
%> $ AmntDaysInfectedToHospital : num 9.5
%> $ AmntDaysNormalToHealthy : num 10
%> $ AmntDaysNormalToIntensive : num 5
%> $ AmntDaysNormalToVentilation : num 3.63
%> $ AmntDaysNormalToDeath : num 5
%> $ AmntDaysIntensiveToAftercare : num 7
%> $ AmntDaysIntensiveToVentilation : num 4
%> $ AmntDaysIntensiveToDeath : num 5
%> $ AmntDaysVentilationToIntensiveAfter : num 30
%> $ AmntDaysVentilationToDeath : num 20
%> $ AmntDaysIntensiveAfterToAftercare : num 3
%> $ AmntDaysIntensiveAfterToDeath : num 4
%> $ GammaShapeParameter : num 1
%> $ FactorPatientsInfectedToHospital : num 0.1
%> $ FactorPatientsHospitalToIntensive : num 0.09
%> $ FactorPatientsHospitalToVentilation : num 0.01
%> $ FactorPatientsNormalToIntensive : num 0.1
%> $ FactorPatientsNormalToVentilation : num 0.001
%> $ FactorPatientsNormalToDeath : num 0.1
%> $ FactorPatientsIntensiveToVentilation : num 0.3
%> $ FactorPatientsIntensiveToDeath : num 0.1
%> $ FactorPatientsVentilationToIntensiveAfter: num 0.7
%> $ FactorPatientsIntensiveAfterToDeath : num 1e-05
%> $ AmntDaysAftercareToHealthy : num 3
%> $ RiskFactorA : num 0.0205
%> $ RiskFactorB : num 0.01
%> $ RiskMale : num 1.5
%> $ AmntDaysIntensiveAfterToHealthy : num 3
%> $ FactorPatientsIntensiveAfterToHealthy : num 0.67
babsim.hospital simulator requires the specification of
arrivalTimesconfpara for the simulation.arrivalTimes, a risk can be specified, i.e., a data.frame with the following entries can be passed to the main simulation function babsimHospital:
time: arrival timeRisk: risk (based on age and gender)Risk values is optional.envs.%> Windows detected. Turning off parallel processing.
First, we illustrate how to generate plots using the simmer.plot package.
In the following graph, the individual lines are all separate replications. The smoothing performed is a cumulative average.
Besides intensiveBed and intensiveBedVentilation, babsim.hospital also provides information about the number of non-ICU beds. The non-ICU beds are labeled as bed.
Summarizing, babsim.hospital generates output for three bed categories:
bedintensiveBedintensiveBedVentilationResources are conceived with queuing systems in mind, and therefore they comprise two internal self-managed parts: aserver, which is the activepart, with a specified capacity that can be seized and released (seeseize); and a priorityqueueof acertain size, in which arrivals may wait for the server to be available
To plot resource usage for three resources side-by-side, we can proceed as follows:
resources <- get_mon_resources(envs)
resources$capacity <- resources$capacity/1e5
plot(resources, metric = "usage", c("bed", "intensiveBed", "intensiveBedVentilation"), items = "server",steps= TRUE)
* Each resource can be plotted separately.
plot(resources, metric = "usage", "bed", items = "server", steps = TRUE)
plot(resources, metric = "usage", "intensiveBed", items = "server", steps = TRUE)
plot(resources, metric = "usage", "intensiveBedVentilation", items = "server", steps = TRUE)
babsim.hospital provides functions for evaluating the quality of the simulation results.babsim.hospital provides a default parameter set, that is based on knowledge from domain experts (doctors, members of COVID-19 crises teams, mathematicians, and many more).fieldEvents <- getRealBeds(data = data$fieldData,
resource = conf$ResourceNames)
res <- getDailyMaxResults(envs = envs, fieldEvents = fieldEvents, conf=conf)
resDefault <- getError(res, conf=conf)
The error is 277.0660283.
babsim plots can be generated.p <- plotDailyMaxResults(res, showBeds = TRUE)
plot(p)
ggplot and plotlycan be used to generate interactive plots.ggplotly(p)
babsim.hospital provides a default parameter set, which can be used for simulations.babsimHospitalPara() provides a convenient way to access the default parameter set:para <- babsimHospitalPara()
#print(para)
para$AmntDaysInfectedToHospital= 5#9.5
para$AmntDaysNormalToHealthy= 7#10
para$AmntDaysNormalToIntensive= 3#5
para$AmntDaysNormalToVentilation= 2#3.63
para$AmntDaysNormalToDeath= 4
para$AmntDaysIntensiveToAftercare= 3#7
para$AmntDaysIntensiveToVentilation= 2#4
para$AmntDaysIntensiveToDeath= 3#5
para$AmntDaysVentilationToIntensiveAfter= 3#6#8#30
para$AmntDaysVentilationToDeath= 4#20
para$AmntDaysIntensiveAfterToAftercare= 3
para$AmntDaysIntensiveAfterToDeath= 4
para$GammaShapeParameter= 1
para$FactorPatientsInfectedToHospital= 0.06#####0.07#0.1
para$FactorPatientsHospitalToIntensive= 0.25#0.09
para$FactorPatientsHospitalToVentilation= 0.005#0.01
para$FactorPatientsNormalToIntensive= 0.3##0.1
para$FactorPatientsNormalToVentilation= 0.003#0.001
para$FactorPatientsNormalToDeath= 0.1
para$FactorPatientsIntensiveToVentilation= 0.1#0.3
para$FactorPatientsIntensiveToDeath= 0.4#0.1
para$FactorPatientsVentilationToIntensiveAfter= 0.76#0.7
para$FactorPatientsIntensiveAfterToDeath= 1e-05
para$AmntDaysAftercareToHealthy= 3
para$RiskFactorA= 0.02048948
para$RiskFactorB= 0.01
para$RiskMale= 1.5
para$AmntDaysIntensiveAfterToHealthy= 3
para$FactorPatientsIntensiveAfterToHealthy= 0.67
babsim provides an interface to optimize the parameter values of the simulation model.conf$simulationDates$StartDate
conf$simulationDates$EndDate
conf$fieldDates$StartDate
conf$fieldDates$EndDate
Sys1 <- Sys.time()
print(Sys1)
library("babsim.hospital")
library("SPOT")
library("simmer")
studyDate <- as.Date( min(conf$simulationDates$EndDate, conf$fieldDates$EndDate ))
resUK <- runoptUK(
expName = paste0("UK-", format(Sys.time(), "%Y-%b.%d-%H.%M-V"), utils::packageVersion("babsim.hospital")),
simData = data$simData,
fieldData = data$fieldData,
TrainSimStartDate = studyDate - 10*7, # "2020-09-02",
TrainFieldStartDate = studyDate - 8*7, # "2020-09-15",
TestSimStartDate = studyDate - 6*7 , #"2020-09-30",
TestFieldStartDate = studyDate - 4*7, #"2020-10-23",
Overlap = 0,
seed = 101170,
repeats = 2,
funEvals = 50,
size = 30,
simrepeats = 2,
parallel = TRUE,
percCores = 0.9,
icu = FALSE,
icuWeights = c(1,2,10),
verbosity = 0,
resourceNames = c("bed", "intensiveBed", "intensiveBedVentilation"),
resourceEval = c("bed", "intensiveBed", "intensiveBedVentilation")
)
Sys2 <- Sys.time()
print(Sys2-Sys1)
runoptUK() returns a list with two elements:
result.df, which is a data.frameresult.df contains the best (optimized) results from the SPOT runs.ukpara as follows:
#ukpara <- resUK[[1]]
#usethis::use_data(ukpara, overwrite = TRUE)
getBestParameter picks out the best and maps it to the variables that are used by the BaBSim.Hospital simulator.print(ukpara)
runopt() optimization from above can be used as follows:para <- getBestParameter(resUK[[1]])
ukpara and can be converted into a babsimhospitalparameter set as follows:para <- getBestParameter(ukpara)
str(para)
%> List of 29
%> $ AmntDaysInfectedToHospital : num 6.09
%> $ AmntDaysNormalToHealthy : num 7.1
%> $ AmntDaysNormalToIntensive : num 6.38
%> $ AmntDaysNormalToVentilation : num 4.89
%> $ AmntDaysNormalToDeath : num 3.21
%> $ AmntDaysIntensiveToAftercare : num 8.67
%> $ AmntDaysIntensiveToVentilation : num 3.6
%> $ AmntDaysIntensiveToDeath : num 6.69
%> $ AmntDaysVentilationToIntensiveAfter : num 27.4
%> $ AmntDaysVentilationToDeath : num 19
%> $ AmntDaysIntensiveAfterToAftercare : num 2.03
%> $ AmntDaysIntensiveAfterToDeath : num 1.3
%> $ GammaShapeParameter : num 0.59
%> $ FactorPatientsInfectedToHospital : num 0.0741
%> $ FactorPatientsHospitalToIntensive : num 0.106
%> $ FactorPatientsHospitalToVentilation : num 0.00625
%> $ FactorPatientsNormalToIntensive : num 0.12
%> $ FactorPatientsNormalToVentilation : num 0.00179
%> $ FactorPatientsNormalToDeath : num 0.0987
%> $ FactorPatientsIntensiveToVentilation : num 0.271
%> $ FactorPatientsIntensiveToDeath : num 0.107
%> $ FactorPatientsVentilationToIntensiveAfter: num 0.798
%> $ FactorPatientsIntensiveAfterToDeath : num 0.0064
%> $ AmntDaysAftercareToHealthy : num 3.48
%> $ RiskFactorA : num 0.484
%> $ RiskFactorB : num 0.0603
%> $ RiskMale : num 1.61
%> $ AmntDaysIntensiveAfterToHealthy : num 4
%> $ FactorPatientsIntensiveAfterToHealthy : num 0.691
conf$simRepeats = 10
res <- modelResultHospital(para=para,
conf=conf,
data = data)
%> Windows detected. Turning off parallel processing.
resOpt <- getError(res, conf=conf)
p <- plotDailyMaxResults(res, showBeds = TRUE)
print(p)
ggplot and plotlycan be used to generate interactive plots.ggplotly(p)